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ABSTRACT 


This  paper  gives  an  algorithm  for  L-shaped  linear  programs  which  arise 
naturally  in  optimal  control  problems  with  state  constraints  and 
stochastic  linear  programs  (which  can  be  represented  in  this  form 
with  an  infinite  number  of  linear  constraints).  The  first  section 
describes  a  cutting  hyperplane  algorithm  which  is  shown  to  be 
equivalent  to  a  partial  decomposition  algorithm  of  the  dual  program. 
The  two  last  sections  are  devoted  to  applications  of  the  cutting 
hyperplane  algorithm  to  a  linear  optimal  control  problem  and 
stochastic  programming  problems. 


L-SHAPED  LINEAR  PROGRAMS  WITH  APPLICATIONS  TO 
OPTIMAL  CONTROL  AND  STOCHASTIC  PROGRAMMING 


by 

R.  M.  Van  Slyke  and  Roger  J.~B.  Wets 

I.  Introduction 

It  has  been  observed  by  many  authors,  see  e.g.  Barr  [2],  Gilbert  [12],  Rosen 
[21,22],  Neustadt  [18],  Whalen  [29,30],  Zadeh  [31],  Pshenichniy  [19],  that  the 
techniques  of  mathematical  programming  can  be  utilized  to  solve  optimal  control 
problems.  The  usual  approach  (although  others  are  possible,  see  e.g.  Dantzlg  [6], 

Van  Slyke  [23])  is  to  discretize  the  system  either  by  finite  difference  approximations 
or  by  considering  the  system  in  sample  data  mode.  If  the  system  dynamics  are  linear 
and  there  are  no  state  space  constraints  various  devices  [6,23]  of  mathematical 
programming  can  be  used  so  that  the  grid  size  or  number  of  sample  points  in  the 
sample  mode  does  not  affect  the  number  of  equations  in  the  associated  mathematical 
program.  This  is  desirable  since  the  computational  effort  for  solving  linear 
programs  by  the  simplex  method  depends  much  more  on  the  number  of  equations 
involved  than  on  the  number  of  variables.  However,  if  state  space  constraints  are 
present  the  number  of  equations  can  grow  astronomically.  This  is  unfortunate, 
especially  in  the  common  situation  where  the  state  space  constraints  are  automatically 
satisfied  for  most  time  periods.  If  confronted  with  problems  of  this  type,  the 
following  heuristic  procedure  suggests  itself:  First  solve  the  problem  without  the 
state  space  constraints,  then  check  if  the  solution  satisfies  all  the  state  space 
constraints.  If  it  violates  some  of  these  constraints,  introduce  only  those-  which 
are  violated  and  solve  this  new  problem.  The  procedure  is  repealed  until  a  feasible 
(and  thus  optimal)  solution  is  attained.  The  algorithm  developed  in  this  paper 
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formalizes  the  ideas  of  this  heuristic  procedure.  Whenever  we  obtain  a  solution 
which  violates  some  state  space  constraint,  we  generate  a  restriction  on  the 
controls  (rather  than  on  the  states)  which  eliminates  the  various  solution  from 
the  feasibility  region. 

The  algorithm  can  be  slightly  modified  to  solve  stochastic  programs  with  re¬ 
course,  first  considered  by  Dantzig  [5]  and  Dantzig  and  Madansky  [7]  under  the  name 
of  two-stage  linear  programs  under  uncertainty.  The  problem  here  is  the  following: 
A  decision  must  be  made  before  the  actual  value  of  some  of  the  parameters  of  the 
problem  is  observed  (it  is  assumed  that  those  parameters  are  known  in  probability). 
Due  to  the  lack  of  knowledge  of  the  particular  outcome  of  the  random  elements  of 
the  problem  discrepancies  may  occur  which,  after  observing  the  actual  values  of 
those  parameters,  are  to  be  corrected  by  selecting  a  particular  recourse  action 
(also  called  second  stage  decision).  One  of  the  difficulties  which  arise  when 
trying  to  solve  such  problems  is  that  a  particular  decision  and  a  particular  out¬ 
come  of  the  random  elements  may  give  rise  to  discrepancies  for  which  there  is  no 
(feasible)  recourse  action.  Thus,  one  should  only  select  decisions,  such  that  for 
every  possible  realization  of  the  random  elements,  a  (feasible)  recourse  action 
can  be  selected  to  correct  the  eventual  discrepancies.  Earlier  treatments  of 
stochastic  programming  ignored  this  difficulty  [7],  [25]  by  assuming  that  the 
structure  of  the  stochastic  program  was  such  that  this  problem  could  not  arise.  It 
did  appear  that  unless  one  made  this  assumption  the  additional  constraints  one  had 
to  introduce  could  be  very  large,  even  infinite  when  the  random  parameters  had 
continuous  distributions. 

In  [7]  Dantzig  and  Madansky  considered  stochastic  programs  with  finitely  dis¬ 
tributed  random  parameters  and  complete  recourse,  i.e.,  for  every  decision  and 
for  every  outcome  of  random  variables  there  exists  a  feasible  recourse.  For 
obvious  practical  reasons  it  seemed  desirable  to  remove  those  restrictive 
assumptions.  The  last  section  of  this  paper  develops  an  algorithm  for  stochastic 
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programs  which  fail  to  satisfy  the  complete  recourse  assumption  as  well  as  the 
finite  distribution  assumption. 

In  [27]  it  was  shown  that  for  stochastic  programs  with  recourse  (random 
right  hand  sides)  the  set  of  feasible  decisions,  represented  by  a  n-vector  x  ,  is 
a  convex  polyhedral  subset  of  Rn  ,  thus  at  most  a  finite  number  of  linear 
constraints  must  be  added  to  the  problem  to  determine  the  set  of  feasible  decisions. 
However,  the  characterization  of  the  feasibility  region  given  in  [27]  is  not  very 
constructive.  The  algorithm  developed  here  generates  these  linear  constraints 
systematically  and  generates  only  those  which  are  violated  by  some  optimal 
decision  candidate,  in  much  the  same  way  as  in  the  control  problem  with  state  space 
constraints . 

The  stochastic  programming  problem  differs  from  the  linear  optimal  control 
problem  in  that  there  is  a  cost  associated  with  the  recourse  actions  which  must  be 
accounted  for.  Dantzig  and  Madansky  [7]  suggest  sampling  to  obtain  the  appropriate 
characteristics  of  the  cost  function  associated  with  the  recourse  problem. 

However,  as  pointed  out  by  Madansky^,  the  utilization  of  sampling  can  lead  to 
inaccuracies.  This  as  we  shall  see,  can  be  avoided  by  using  a  gradient  method  rather 
than  a  cutting  plane  method. 

In  Section  2, an  algorithm  which  is  essentially  the  same  as  the  algorithm 
developed  by  Benders  [3]  is  described  and  a  geometric  interpretation  is  given. 
Section  3  exhibits  the  duality  between  this  algorithm  and  a  variant  of  the 
decomposition  algorithms  of  Dantzig  and  Wolfe.  The  applications  of  this  algorithm 
to  optimal  control  problems  with  state  space  constraints  and  stochastic  programs 
with  recourse  are  developed  in  Sections  4  and  5  respectively.  Now,  let  us  give  a 
mathematical  formulation  of  the  linear  program  we  are  interested  in. 


SIGMAP  Conference  on  Stochastic  Programming,  Princeton,  New  Jersey*  (December  1965). 

4-4- 

This  was  pointed  to  us  by  E.  Balas. 


Thus,  in  order  to  speed  up  computation  and  limit  storage  requirements  it  is 
desirable  to  work  mainly  with  the  constraints  (l.a)  and  consider  the  constraints 
(l.b)  and  the  variables  of  y  only  when  needed. 
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2.  A  Cutting  Plane  Algorithm 
tu _ Feasibility. 

2 

Instead  of  problem  (1),  let  us  first  consider  the  special  case  where  c  ■  0  . 
(This  corresponds  to  the  problem  arising  in  optimal  control  problems  with  state 
space  constraints) : 


(2) 


Minimize  z 
subject  to 


1 

C  X 


A11 

A  x 


-  b 


a21  i  a22  ,2 

A  x  +  A  y  ■  b 


x  >  0,  y  >  0 


In  this  case  the  algorithm  proceeds  as  follows.  First  solve  the  simpler  linear 
program 

Minimize  z  »  c*x 

(3)  subject  to  A^x  ■ 

x  >  0 

whose  optimal  solution  we  denote  by  x  .  For  the  time  being  we  assume  that  (3)  is 
solvable.  If  x  satisfies  the 

22  2  21- 

Feasibility  Criterion :  There  exists  y  >  0  such  that  A  y  »  b  -  A  x  , 

then  x  and  some  y  determine  a  feasible  (and  thus  optimal)  solution  to  (2).  We 
denote  by  the  set  of  all  x  satisfying  the  feasibility  criterion. 

If  x  does  not  satisfy  the  Feasibility  Criterion,  we  generate  a  constraint 
involving  only  x  which  is  violated  by  x  but  satisfied  by  any  feasible  solution  to 
(2).  This  constraint  is  then  added  to  the  constraints  of  problem  (3).  This  added 
constraint  has,  in  a  sense  to  be  made  precise  later,  (Section  2,  E)  the  property 
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that  it  cuts  deepest  into  the  set  =  {x|a^x  *  b\  x  >  0}  .  The  process  is 
then  repeated  until  an  optimal  solution  to  the  augmented  problem  (3)  satisfies 
the  feasibility  criterion.  We  will  show  that  we  have  to  add  at  most  a  finite 
number  of  constraints  to  (3)  in  order  to  achieve  this  goal. 

To  determine  whether  x  satisfies  the  Feasibility  Criterion  or  not  we  try 
to  find  a  nonnegative  solution,  y  ,  to 


(4) 


* 2  2 
A  y 


A 


21- 
x  . 


22  22 

This  can  be  considered  geometrically.  Let  pos  A  ={t|t=A  y  ,  y  >  0}  be 

22 

the  closed  convex  cone  generated  by  the  columns  of  A  .  Then  x  satisfies 

2  21-  2? 

the  feasibility  criterion  if  and  only  if  b  -  A  x  c pos  A  .  If  not,  i.e.,  if 
2  21-  22 

b  -  A  x  t  pos  A  ,  there  is  a  hyperplane  through  the  origin  separating  strictly 
2  21-  22 

b  -  A  x  and  pos  A  .  Such  a  hyperplane,  say  {x|ox  =  0}  ,  is  determined  by  its 

22  2  21  - 

normal  o  which  must  satisfy  at  <  0  for  t  e  pos  A  and  o[b  -  A  ]x  >  0  . 


Figure  1. 

The  normals  a  ,  which  are  needed,  are  generated  using  a  slight  variant  of 


the  Phase  I  procedure  for  the  simplex  method.  We  solve 
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(5) 


Minimize  w  =  ev+  +  ev 

22  +  _  2  21- 
subject  to  A  y+Iv  -Iv  *  b  -  A  x 

y  >  0  ,  v+  >  0  ,  v  >0, 


where  e  is  a  row  vector  of  l's,  I  is  a  x  m2  identify  matrix  and  v  , 
and  v  are  n^-vectors  of  variables. 

Problem  (5)  has  always  an  optimal  solution  with  w  >  0  .  x  satisfies  the 
feasibility  criterion  if  and  only  if  at  the  optimal  w  •  0  .  If  at  the  optimum 
w  >  0  ,  then  there  exist  dual  variables  0  satisfying 


<  0 


(6)  -e  <  0  <  e 

0  [b^  -  A^x]  -  Minimum  w  >  0  . 

Thus  o  has  the  desired  properties.  In  the  next  sections  we  show  that  the  o's 
generated  by  solving  (5)  are  optimal  in  some  sense  and  give  the  geometrical 
interpretation  in  more  detail. 

2  21 

In  order  for  x  to  be  feasible  it  is  clear  that  b  -  A  x  must  be  on  the 

22 

same  side  of  the  hyperplane  {t|ot  *=  0}  as  pos  A 
Thus,  x  feasible  implies  that 

2  21 

o [b  -  A  x]  <  0  . 

Thus,  we  add  the  constraint 

(7)  [oA^]  x  >  ob^ 
to  the  linear  program  (3). 

It  is  also  possible  that  when  solving  problem  (3)  (or  even  after  a  few 


additional  constraints  have  been  added)  we  discover  that  (3)  is  unbounded.  Thus, 
the  solution  to  (3)  is  no  longer  given  in  terms  of  a  particular  vector  x  but 
we  are  given  a  half-line  in  ,  say  x  +  Xxc  ,  X  >  0  ,  on  which  cx  decreases 


monotonically  to  -  00 

as 

X  goes  to  +  “>  . 

We 

have  that: 

(8)  Proposition 

Tf  A21* 

If  -A  x 

c 

and 

2  21- 

b  -  A  Xp  belong 

to 

22 

pos  A  ,  then 

(2) 

is 

unbounded. 

If 

A21"  J  A22 

-Ax?  pos  A 
c 

then 

every  solution 

to 

(2) 

must  satisfy  the  constraint. 


(9)  [oA^]x  >  ob^ 

which  is  violated  by  xp  +  Xx^  for  X  sufficiently  large,  where  o  denot 
the  vector  of  optimal  simplex  multipliers  corresponding  to  the  optimal 
solution  to 

Minimize  w  =  ev+  +  ev 


22  +  -  21- 
subject  to  A  y  +  Iv  -  Iv  =  -  A  x 


y  >  0  ,  v  >0,v  >0. 

21-  22  2  91-  22 
If  -A  x^  e  pos  A  but  b  -A"  x^  i  pos  A  then  every  feasible 

solution  to  (2)  must  satisfy  the  constraint  generated  by  solving  the 

linear  program  (5)  where  x  is  set  equal  to  x^  . 


Proof : 


21-  2  21-  22 

The  conclusion  is  immediate  if  -A  x  and  b  -  A  x  belong  to  pos  A 

c  P 

21-  22  -  2  21-  21-  22 
If  -A  Xq  i  pos  A  then  for  some  X  ,  b  -  A  x^  -  XA  xc  ^  pos  A  for  all 

X  >  X  .  To  see  this,  it  suffices  to  observe  that  if  u  are  the  optimal  simplex 

22  21- 

multipliers  for  (10)  then  oA  •  0  and  oA  x  -  0  .  Thus,  bv  scle  :* in>  X 

2  21  -  21- 

sufficiently  large  o (b  -  A  x  -  XA  x^ )  can  also  he  made  arbitrarilv  small. 
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-  2  21-  - 
Set  X  =  0  if  o(b  -A  Xp)  <  0  ,  otherwise  select  X  such  that 

2  21-  -  21- 

o(b  -  A  x  -XA  x)“0.  Then  for  all  X  >  X  ,  o  determines  a  hyperplane 
22  2  21-  21- 

separating  pos  A  and  (b  -  A  Xp  -  XA  xc>.  It  follows  that  every  x  in 


K,  such  that  x  *  (x  +  Xx  )  +  px 
1  p  c  c 


p  >  0  violates  (9)  which  must  be  satisfied 


21-  22  2  21-  22 

by  every  feasible  solution  to  (2).  If  -A  xc  c  pos  A  but  b  -  A  Xp  ^  P°8  ^ 

2  21-  21-  22 
either  b  -  A  x  -  XA  xc  does  not  belong  to  pos  A  for  all  X  >  0  or  there 

-  2  21-  21-  22 
exists  X  such  that  if  X  >  X  ,  b  -  A  x^  -  XA  x£  belongs  to  pos  A  .  Now 

let  o  denote  the  optimal  simplex  multipliers  obtained  from  (5)  by  setting 

2  21-  21-  22 
x  *  Xp  .  If  for  all  X  >  0  ,  b  -  A  x^  -  XA  xc  does  not  belong  to  pos  A  , 

the  ray  Xp  +  Xx£  violates  for  all  X  the  constraint  (7)  so  generated  and  thus 

this  particular  extreme  ray  is  eliminated  from  the  feasible  solution.  On  the 
2  21-  21-  22 

other  hand,  if  b  -  A  Xp  -  XA  x£  belongs  pos  A  for  A  >  A  ,  the  points 


2  21  -  —  21- 

b  -  A  (x  +  Xx  )  -  pA  x  p  >  0 

p  c  c 

satisfy  the  constraints  and  the  ray  (Xp  +  Xx^)  +  px^  has  not  been  eliminated  from 
the  set  of  feasible  solutions. 

We  can  thus  summarize  the  procedure  to  find  an  optimal  solution  to  (2) ,  as 
follows  : 

If  (3)  (with  as  without  additional  constraints)  is  solvable  with  x  *  x  .  We 
then  solve  (5).  If  w  <■  0  then  x  is  an  optimal  solution  for  (I).  Otherwise  we 
generate  a  constraint  of  the  type  (7)  which  is  then  added  to  the  constraints  of  (3). 
If  (3)  (with  or  without  additional  constraints)  is  unbounded  with  a  direction  of 
decrease  for  cx  given  by  x  ■  Xp  +  Xxc  ,  X  >  0  .  We  then  solve  (10),  and  (5) 
with  x  ■  Xp  .  Let  w  and  w  denote  the  optimal  value  for  (10)  and  (5) 
respectively.  If  w  ■  w  -  0  then  (2)  is  unbounded.  If  w  >  0  ,  we  use  the 
optimal  multipliers  of  (10)  to  generate  a  constraint  of  the  type  (9)  which  is 
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added  to  the  constraints  of  (3);  if  w  =  0  but  w  >  0  ,  we  generate  a  constraint 
of  the  type  (7). 

Clearly  this  process  is  finite  since  each  o  corresponds  to  a  basis  for  (5) 
(or  (10))  of  which  there  is  a  finite  number  and,  moreover,  no  constraint  can  be 
repeated.  Obviously,  no  constraint  of  type  (9)  will  be  generated  after  we  obtain 
bounded  solution  to  (3). 

It  is  conceivable,  of  course,  that  the  number  of  bases  of  (5)  or  (10), 
corresponding  to  a  particular  a  could  be  very  large,  so  chat  the  number  of 
generated  constraints  could  be  large  compared  to  the  number  of  original  constraints 
(l.b)  in  which  case  the  proposed  algorithm  might  be  inefficient.  However,  since  we 
only  add  binding  constraints  which  have  a  deepest  cut  property  (as  we  shall  see 
later)  and  if  properties  (i),  (ii)  and  (iii)  mentioned  in  the  introduction  are 
satisfied,  this  seems  unlikely. 

Another  useful  property  of  the  algorithm  is  that  in  adding  new  constraints  to 
(3),  the  next  iteration  already  has  a  basic  solution  which  is  infeasible  only  for 
one  basic  variable.  The  basis  is  the  basis  for  the  previous  iteration,  plus  the 
slack  variable  for  the  added  constraint.  Thus,  each  successive  x  can  be  easily 
obtained  by  a  few  steps  of  the  dual  simplex  method. 

B.  Optimality 

2 

We  now  return  to  our  original  problem  (1),  i.e.,  to  the  case  when  c  may 
be  different  of  zero.  Obviously,  problem  (1)  is  equivalent  to: 

Minimize  c^x  +  0 

(11)  subject  to  Q(x)  <  0 

x  t  K  =  K]  ^  K2 

where 

(12)  Q(x)  =  {Min  c^ylA^y  =  b^  -  A^x  ,  x  >  0)  . 


J 
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We  first  observe  that 

(13)  Proposition: 

For  all  x  e  ,  Q(x)  is  either  a  finite  convex  function  or 
Q(x)  is  identically  -  00  . 

Proof : 

For  all  x  e  ,  the  linear  program 

2 

Minimize  c  y 

22  2  21 

(14)  subject  to  A  y  -  b  -Ax 

y  >  0 

is  feasible.  Moreover,  for  all  x  c  ^  .  (14)  is  unbounded  if  and  only  if  the 
22 

linear  system  ttA  <  c  is  inconsistent.  Thus,  if  (14)  is  unbounded  for  some  x  , 
it  will  be  unbounded  for  all  x  .  It  remains  to  show  that  if  Q(x)  is  finite  on 
the  convex  set  ,  then  it  is  convex.  Consider  x^,  x1  e  ^  and 
x*  *  (l-A)x^  +  Ax1  where  Ae[0,l]  ,  and  let  y^,  y1  and  yA  be  optimal  solution 
to  (14)  when  x  equals  x^,  x1  or  xA  respectively.  Then, 

(l-A)Q(x°)  +  AQ(x1)  +c[(l-A)y°  +  Ay1]  >  cyA  +  Q(xX) 

since  (l-A)y^  +  Ay1  is  a  feasible  solution  (but  not  necessarily  optimal)  to  (14) 
when  x  equals  xA  . 

Moreover, 

(15)  Proposition: 

Suppose  Q(x)  is  finite,  let  n  denote  the  optimal  simplex 
multipliers  corresponding  to  the  solution  of  (14)  with  x  =  x  , 


then  the  linear  function 
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(16)  (ttA^)x  -  (irb^) 

Is  a  support  of  Q(x)  . 

Proof : 

Since  7t  is  optimal  for  (14)  with  x  =  x  ,  then  by  duality  theory  for  linear 
programming  we  have  that 


tt  (b ^  -  A^x)  =  Q(x)  . 

By  assumption  Q(x)  is  finite  and  thus  for  all  x  e  K,  ,  j  is  a  feasible  solution 
for  all  duals  of  (14): 

(17)  Maximize  tt  (b^  -  A^x) 

22  2 
ttA  <  c 

but  it  is  not  necessarily  an  optimal  solution.  Thus,  again  by  duality  theory  we 
have 


7T  (b  ^  -  A^x)  <  {Max  tt  (b  ^  -  A^x)|nA^  <  c^}  =  Q(x)  for  all  x  t  . 

Even  though  the  following  observation  is  not  absolutely  necessary  for  the 
subsequent  development,  it  is  worthwhile  to  note  that 

(18)  Proposition: 

Suppose  Q(x)  is  finite  on  K7  ,  then  Q(x)  is  a  convex 
polyhedral  function. 

Proof  : 

By  letting  x  range  over  ,  we  see  that  only  a  finite  number  of  supports 

to  Q(x)  of  the  type  (16)  can  be  generated,  since  every  v  corresponds  to  a 
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particular  basis  of  A  and  A  has  only  a  finite  number  of  square  nonsingular  sub- 
matrices.  Moreover,  for  all  x  c  there  is  some  support  of  type  (16)  which 

meets  Q(x)  at  x  .  Thus,  the  upper  envelope  of  this  finite  number  of  linear 
supports  coincide  with  Q(x)  . 

The  process  to  obtain  an  optimal  solution  to  (1)  or  equivalent  to  (11)  is 
very  similar  to  the  one  already  described  for  finding  a  feasible  solution.  Suppose 
x  is  a  feasible  solution,  i.e.,  x  e  K  ■  ,  and  (14)  is  solvable  with 

x  -  x  .  Let  tt  be  the  corresponding  optimal  simplex  multipliers.  Then: 

Q(x)  ■  n(b2  -  A2^x)  . 

Moreover,  by  convexity  of  Q(x)  and  the  properties  of  tt  given  in  (15),  it  follows 
that 

Q(x)  >  irb2  -  [wA21]x 

for  all  x  in  K  .  Thus,  a  pair  (x,0)  is  feasible  for  (11)  only  if 

-  2  -  21 
8  >  irb  -  [  nA  ]x 

which  we  can  also  write 


(19)  [wA2^]x  +  0  >  Tib2 

On  the  other  hand  if  (x*\  8^)  are  optimal  for  (11)  and  are  the  optimal 

simplex  multipliers  obtained  from  (14)  by  substituting  x  for  x^  ,  we  have  that 


Q(x°) 


02  0.21  0 
TT  b  -  IT  A  X 


The  optimality  of  x^  implies  that  cx  +  Q(x)  cx°  +  Q(x^)  for  all  x  In  K  . 
From  0^*  >  Q(x^)  and  0  unrestricted  in  (11)  it  follows  that  0^  =  Q(x^)  . 
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These  two  last  observations  allow  us  to  construct  a  finite  procedure  for 
finding  an  optimal  solution  to  (1).  Say  (x  ,  0  )  is  an  optimal  solution  to  the 
linear  program 

(20)  Minimize  c^x  +  0 

(20. a)  subject  to  [ti^A^Ix  +  0  >  (n *"b^ )  H*l,  ...»  k-1 

x  c  Kr  ^  K2  . 
k 

We  then  solve  (14)  with  x  =  x  .  If  (14)  is  unbounded  then  (1)  is  unbounded.  If 
k+1 

not,  let  tt  denote  the  optimal  simplex  multipliers.  Then: 

(21)  Optimality  Criterion: 

If  k+l  ?  ?  1  k  k 

If  e  =  tt  [t  -  Ax]  then  x  is  an  optimal  solution  to  (1).  If 

k  k 

0  <  Q(x  )  we  add  the  constraint 


r  k+1  21,  A  k+1.  2 

[tt  a  ]x  +  0  >  tt  b 


to  the  constraints  of  (20),  which  has  the  effect  of  eliminating  the  solution 
k  It 

(x  ,  0  )  from  the  set  of  feasible  solutions  of  (2).  The  algorithm  is  initiated 

with  x^  minimizing  c^x  on  K  and  0*^  =  -  u  . 

Now  suppose  that  (20)  is  unbounded  after  at  least  one  constraint  of  type  (20. a) 

has  been  introduced.  Note  that  in  such  a  case,  (20)  cannot  be  unbounded  for  some 

fixed  x  and  0  =  -  00  ,  since  0  must  satisfy  the  constraint  of  type  (20. a).  Thus, 

there  exists  some  ray,  say  x  +  Ax  ,  A  '•  0  on  which  the  oblective  of  (20)  can  >e 

p  c 

22 

pushed  to  -  °°  .  Checking  if  this  ray  belongs  to  pos  A  has  been  dealt  with  in 
the  previous  section.  If  not,  we  generate  constraints  of  type  (7)  or  (9).  Now 


2  21  21  22 
suppose  b  -  A  Xp  and  -AA  x^  belong  to  pos  A 

solution  to  the  linear  program: 


Let  y^  be  an  optimal 
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Minimize  c  y 

22  21 

subject  to  A  y  ■  -  A  xc 
y  >  0 

and  let  tt  be  the  corresponding  vector  of  optimal  simplex  multipliers.  If 

12  12 
c  x  +  c  y  <0  then  obviously  (1)  is  unbounded.  If  cx  +  cy  >0  then 
c  c  c  c 

xc  is  not  a  desirable  unbounded  direction  since  letting  X  go  to  +  00  in 
Xp  +  Xxc  would  push  the  objective  of  (1)  to  +  00  .  In  this  case  adding  the 
constraint 

[irA^]x  +  0  >  [irb^] 

to  (20)  would  eliminate  the  direction  xc  from  the  desirable  (optimal)  solutions 
1  2 

of  (20).  If  c  x£  +  c  yc  ■  0  then  no  point  of  the  ray  x^  +  Xx£  will  be 

preferable  to  x^  as  a  solution  to  (1),  thus  adding  the  above  constraint  to  (20) 

will  keep  x  in  the  set  of  feasible  solutions  of  (20)  but  will  eliminate  the 
P 

other  points  of  the  ray. 

22 

This  process  is  obviously  finite  since  each  n  corresponds  to  a  basis  of  A 
and  these  are  finite  in  number  Moreover,  no  n  can  be  generated  twice  since  this 
would  lead  to  a  constraint  already  present  which  could  not  be  violated  by  the  solution 
at  hand.  In  this  section  we  have  assumed  that  each  x  generated  is  a  feasible 
solution,  if  x  i  then  one  may  have  to  introduce  constraints  of  type  (7)  or  (9) 

before  continuing  the  search  for  an  optimal  solution  to  (1) . 

C.  Summary  of  the  Algorithm. 

Step  1: 

Solve  the  linear  program 


i 
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(22) 

Minimize  z  =  c^x  +  0 

(22. a) 

subject  to  A^x  =  b^ 

(22. b) 

[okA21]x  >  [A2] 

k=l ,  ....  s 

(22. c) 

[nkA21]x  +  0  >  hkb2] 

k=l ,  ....  t 

x  >  0  . 

Initially,  s  =  t  =  0  .  6  is  set  equal  to  -  00  and  is  deleted  from  the 

actual  computations  as  long  as  there  are  no  constraints  of  type  (22. c).  If  (22)  is 
infeasible  so  is  (1)  and  we  terminate.  If  (22)  is  solvable  go  to  Step  2  if  (22)  is 
feasible  but  unbounded  go  to  Step  2*. 

Step  2: 

0  £ 

(22)  is  solvable.  Let  (x  ,  6  )  be  an  optimal  solution  to  (22).  Use  the 
simplex  method  (Phase  I,  Phase  II)  to  solve. 


(23) 


Minimize  w  = 

22 

subject  to  A  y  = 

y  > 


2 

c  y 


.2  a21  * 

b  -  A  x 


0  . 


If  (23)  is  infeasible,  i.e..  Phase  I  terminates  with  the  infeasibility  form 
different  of  zero,  we  use  the  multipliers  so  generated  to  construct  a  constraint  of 

the  form  (22. b).  If  (23)  is  feasible  and  unbounded  so  is  (1)  and  we  terminate.  If 

£  £  £ 

(23)  is  solvable  and  Min  w(x  )  =  6  then  is  x  is  optimal  and  we  terminate. 
Otherwise,  we  use  the  multipliers  so  generated  to  construct  a  constraint  of  the  form 
(22. c)  and  return  to  Step  1. 


Step  2 1 : 


(22)  is  feasible  but  unbounded. 


Let 


X  >  0  be  a  ray  of  unbounded 
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decrease  of  c^x  .  We  then  solve  (23)  with  -  A^x^  replaced  by  -A^x^  . 

If  this  problem  is  infeasible  (i.e.f  Phase  I  terminates  with  positive  objective 

value),  we  use  the  optimal  simplex  multipliers  to  generate  a  constraint  of  type 

l  £ 

(22. b).  If  this  problem  is  feasible  let  y^  be  the  optimal  solution  and  tt  the 

i.  £ 

associated  simplex  multipliers.  Now  solve  (23)  with  x  *  x  If  this  new 

P 

problem  is  infeasible,  we  generate  a  constraint  (22. b)  as  in  Step  2.  Otherwise, 

1  £  2  5 

either  cx  +  c  y  '  <  0  in  which  case  (1)  is  unbounded  and  we  terminate  or 
c  c 

1  £  2  £  £ 

c  x£  +  c  yc  >  0  and  then  tt  is  used  to  generate  a  constraint  of  type  (22. c) 
we  return  to  Step  1. 

Finally,  it  is  not  difficult  to  see  that  if  so  desired  (e.g.  in  order  to  keep 

the  data  related  to  problem  (22)  in  the  easy  access  memory)  it  is  possible  to  remove 

those  constraints  of  (22. b)  and  (22. c)  which  are  slack,  although  they  may  be 

generated  again  and  have  to  be  re-introduced.  This  also  necessitates  a  new 

finiteness  proof  which  is  based  on  the  fact  that  upon  taking  suitable  account  for 

degeneracy,  the  objective  value  c^x  +  Q(x)  corresponding  to  every  feasible 

solution  to  (23)  -  generated  in  Step  2  -  is  monotonically  decreasing  so  there  are 

only  a  finite  number  in  which  (23)  has  a  feasible  solution.  On  the  other  hand, 

between  feasible  solutions  to  (23)  when  constraints  of  the  form  (22. b)  are  being 

1  k 

introduced,  the  value  of  c  x  is  monotonically  increasing  so  that  a  feasible 
solution  to  (23)  always  occurs  after  a  finite  number  of  steps. 

D.  Some  Geometric  Characterizations 

We  have  already  pointed  out  that  checking  if  a  particular  point,  say  x  ,  is 

2  21- 

feasible  corresponds  to  determining  if  b  -.A  x  belongs  to  the  cone 
22  22 

pos  A  ■  {t|t  ■  A  y  ,  y  >  0}  .  Similarly,  if  at  some  stage  the  program  (22) 

2  21  e 

yields  an  unbounded  direction,  then  solving  (23)  with  b  -  A  x  replaced  with 

21  £  21  £ 

-A  xc  corresponds  to  determining  if  the  ray  \(-A  xc>  ,  X  >  0  belongs  to  the 
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22  l  l 

cone  pos  A  .  Even  checking  for  optimality  of  a  given  pair  (6  ,  x  )  can  be 


/  6l  \ 

viewed  as  determining  if  (  ^2  ^21^4,  /  belongs  or  does  not  belong  to  the  cone 

pos  (>)■  |  t  =  c^y  ,  A^y  ,  y  >  o|  •  In  this  section  and  the  following  one, 

we  limit  our  discussion  to  the  case  when  checking  for  feasibility,  i.e.,  x  in  , 
but  in  view  of  the  above  observations  our  remarks  can  be  adapted  equally  well  to 
the  other  parts  of  the  algorithm. 

Suppose  x  e  ,  but  does  not  belong  to  ^  •  Then  solving  the  linear 

program  (5)  yields  w  >  0  .  At  the  same  time  we  generate  some  o  ,  which  corresponds 

to  a  particular  basis  of  the  matrix  (A^,  I,  -I)  .  The  basic  solution  contains  at 

least  one  artificial  variable,  i.e.,  a  component  of  the  vector  (v+,  v  )  ,  at 

2  21  22 

positive  level.  Since  otherwise  w  =  0  and  b  -Axe  pos  A  .  We  have  that 


Proposition : 

Suppose  the  optimal  solution  to  (5),  with  w  >  0  ,  contains  exactly 

one  artificial  variable.  Then,  a  is  the  normal  of  a  supporting 

22  22 
hyperplane  of  pos  A  determining  a  (n^-^-dimensional  face  of  pos  A 


Proof : 


22 

First  note  that  this  (n^-lJ-dimensional  face  may  be  pos  A  itself  viz. 

22 

if  pos  A  is  of  dimension  (n^-l).  Also,  by  the  hypothesis  of  this  proposition 
22 

the  cone  pos  A  has  at  least  dimension  n^-l  •  By  assumption,  there  are  at 

22  22  22  th 

least  n^-l  columns  of  A  such  that  oA^  =  0  where  A^  denotes  the  j 

22 

column  of  A  .Of  these,  n^-l  are  linearly  independent  since  ir^-l  belong  to 

9  9  9  9 

the  basis.  Let  F  =  { t  |  t  =  T  A.,yj  ,  y  >  0}  where  J  =  {j  |  oA  =  0}  .  It 

JW  J  22  ,  *1 

now  suffices  to  observe  that  F  =  pos  A  ^  (t  j  ot  =  0)  ,  that  {t  |  o  =0}  is 

22 

a  supporting  hyperplane  of  pos  A  and  that  F  has  dimension  n^-l  since  it 
contains  n^-l  linearly  independent  points. 

Thus,  if  it  is  possible  to  obtain  a  solution  to  (5)  with  only  one  artificial 
variable  in  the  optimal  basis  it  follows  that  o  determines  a  (n^-l)  dimensional 
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22 

faces  of  pos  A  to  construct  the  constraints  (22. b)  one  expects  that  fewer 

22 

need  be  generated.  In  particular,  the  number  of  deficiency  1  faces  of  pos  A 

22 

Is  much  smaller  than  the  number  of  basis  of  (A  ,  I  ,  -I)  In).  However,  It  Is 

22 

not  always  possible  to  obtain  a  (n^-l)  face  of  pos  A  .  In  fact  as  is 
Indicated  in  the  next  proposition,  it  is  sometimes  possible  to  obtain  solutions  to 
(5)  such  that 


{t  |  ot  ■  0}  n  pos  A^  ■  {0}  . 

(25)  Proposition; 

2  21 

Suppose  x  i  K2  and  (b^  -  A^x)  is  different  of  zero  for  all 
1  ,  and  that  for  all  j 


(26) 


l 


{i|(bj  -  A^x)  >  0)  -  (i|  (bj  -  A*Jx)  <  0} 


.  22 
AiJ  > 


l 


21 


22 

ij 


22  22 

holds,  then  no  column  A^  of  A  will  figure  in  the  optimal  basis  of 
(5).  A^  denotes  the  iC^  row  of  A^  . 

Moreover,  one  should  realize  that  the  conclusion  of  the  above  proposition 
depends  very  much  on  the  selection  of  cost  coefficients  +1  for  the  artificial 
variables  in  the  lnfeaslbl llty  form.  In  fact,  any  set  of  positive  numbers  could 
be  selected  as  cost  coefficients  for  the  infeasibility  form.  Thus,  an  obvious 
complement  to  (25)  is 


(27)  Corollary; 

2  21 

Suppose  x  i  Kj  and  (b  -  A^x)  is  different  of  zero  for  all  i  , 

and  for  all  j  and  all  sets  of  positive  numbers  w,  ,  . . . ,  u 

1  m2 


the  relation 


20 


(28) 


{ i  I  (b 


Ajix)  >  Ot 


.22 

uiAij  > 


21 


,22 

ulAlj 


{ i  I  (b±  -  a“x)  <  0} 


22  22 

holds,  then  no  column  of  A  will  figure  in  the  optimal  basis 

of  (5).  A^  denotes  the  i^  row  of  A^  . 

To  see  that  the  condition  (28)  is  not  vacuous  ,  consider 

A22-(22)  and  b2-  A21,  -(:}). 

Obviously  the  condition  (28)  is  much  weaker  than  (26)  since  it  allows  for 

some  perturbation  of  the  coefficients  of  the  objective  function  in  (5).  It  also 

indicates  how  one  may  modify  (5)  in  order  to  be  able  to  increase  the  number  of  tlu 
22 

columns  of  A  figuring  in  the  optimal  basis.  This  would  mutually  increase  the 

22 

dimension  of  the  face  of  pos  A  determined  by  the  corresponding  o  .  In 
practice,  this  would  involve  a  parametric  study  of  the  linear  program  (5).  The 
constraints  (22. b)  so  generated  would  generally  be  "better"  then  those  obtained 
by  solving  (5)  but  whether  the  extra  computation  is  justified  can  probably  be 
discovered  only  by  experience  in  using  the  algorithm. 

E.  A  "Deepest  Cut"  Property 

As  we  mentioned  earlier  the  constraints  (7)  obtained  by  solving  (5)  have  a 

deepest  cut  property  with  interesting  geometrical  interpretations  which  we  now 

examine.  The  linear  program  (5)  can  be  interpreted  as  finding  the  nearest  point 
22  2  21 

in  pos  A  tod=b-Axin  the  sense  of  the  2^  norm,  i.e.. 


Min 


z  -  d 


!  I 


(29) 


22  m2 

subject  to  Zb  pos  A  R 


2  21- 
d  =  b  -  A  x 


wher 


j  denote  the  2  norm  given  by 


m_ 

2 1  i  ■  .Vzi 
1  =  1 


21 


«n„ 


The  £.^  norm  is  defined  on  the  space  R  of  column  n^-vectors .  Associated 
m2  m2  * 

with  R  is  its  dual  space  (R  )  which  may  be  identified  with  all  real  valued 


m,. 


linear  functions  on  R  .As  is  well  known,  any  linear  function  f(z)  in 


can  be  represented  in  a  one  to  one  way  as  a  matrix  product  ir.x  of  a  m2 


dimensional  row  vector  tt  ,  and  the  column  vector  x  .  We  shall  thus  think  of 

m2  *  „m2 
(R  )  as  a  space  of  row  vectors  with  the  same  dimension  as  R  .A  hyperplane, 

m2 

H  ,  passing  through  the  origin  of  R  can  be  represented  in  the  form 

m2  * 

H  *  { z  |  -rr  z  -  0}  for  some  tt  i  0  in  (R  )  .  However,  this  representation  is  not 
unique  since  (6n) z  =  0  determines  the  same  hyperplane  for  any  real  number 


6^0.  To  resolve  this  ambiguity,  we  specify  that 


-  1  .  Where 


m 


2  * 


is  a  norm  defined  on  (R  )  This  norm  can  be  defined  quite  naturally 


m. 


[9],  using 


on  R  by  means  of  the  following  relation 


(30) 


|  tt  ||-  max  { 7tz 


z 


'  1}  . 

m. 


It  is  easily  seen  that  | |  n  | | A  -  | |  n  | |^  -  max  |  |  ,  the  norm  .  A 

given  n  ,  determines  also  a  half  space  S  =  {z  |  ttz  <  0}  which  is  bounded  by  H  . 
The  condition  that  |  |  rr  |  |^  <=  1  and  the  specification  of  which  half  space  is  to  be 
determined  uniquely  defines  tt  . 

The  dual  of  (5)  is 


(31) 


Max  o  (b^  -  A^x] 

22 

subject  to  oA  <0 

k1!  <  1  i  -  1,  .. .  ,  m 


Let  us  interpret  (31)  in  the  language  developed  here.  <1  Is  a  m2 

m2  * 

dimensional  row  vector  which  is  an  element  of  (R  )  It  determines  a  half 

2  2 

space,  S  ,  by  S  =  lx  |  ox  <  0}  which  includes  pos  A  .  This  follows  from 

22  22 

the  relation  z  e  pos  A  implies  that  z  =  A  y  for  some  y  >  0  ;  hence 


22 


22 

z  *  (a A  )y  <  0  .  The  relations  |o^|  <1  i  ■  1 .  is  equivalent  to 


ii  0m2  * 

| |  a  | |  <  1  .  Of  all  elements  of  (R  )  satisfying  these  conditions  we  are 

2  21- 

to  find  one  which  maximizes  o[b  -  A  x]  .  Let  us  now  examine  the  geometrical 

2  21- 

interpretation  of  maximizing  o[b  -  A  x]  . 

The  distance  from  a  point  z  to  a  hyperplane  H  given  by  H  =  {z|oz  *  01 
or  equivalently  from  the  origin  to  the  plane  H-  =  {z|oz  =  oz)  can  be  obtained 
by  solving  the  linear  program: 


Minimize  ez  +  ez 


subject  to  oz  -  oz 


oz 


z  ,  z  >  0  , 


The  optimal  solution  is  obviously  determined  by  »  z^ 


o  z  ; 


z1  -  z\  -  °* 


z  to  H  is 


for  i  4  v  where 


max  | o , 

i 


Thus,  the  distance  from 


max  o  . 
l 


oz  = 


1 

i  - 1 

TT 

|  oz 

IL 

5)) 

can 

Thus,  problem  (31)  (i.e.,  the  dual  of  (3))  can  be  interpreted  as  finding 

m2  *  22 
oc(R  )  determining  a  supporting  hyperplane  of  pos  A  which  is  as  far  as 

2  21- 

possible  from  b  -  A  x  in  the  sense  of  the  norm.  Moreover,  by  the  duality 

theory  of  linear  programming,  we  have  that  this  maximum  distance  is  equal  to  the 
2  21-  22 

>1  distance  of  b  -  A  x  from  pos  A  .  Thus,  in  terms  of  the  H  norm  we 
have  generated  a  "deepest  cut." 
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3.  The  Partial  Decomposition  Algorithm 

A  very  natural  approach  to  L-shaped  programs  Is  via  the  decomposition 
algorithm  of  Dantzig  and  Wolfe  [8].  Nonetheless,  if  (1)  has  the  properties 
mentioned  in  the  introduction,  the  straightforward  application  of  the 
decomposition  algorithm  to  problem  (1)  does  not  take  advantage  of  the  structure 
of  the  problem. 

Decomposition  can,  however,  be  advantageously  applied  to  the  dual  of  problem 


(1): 


1  2 

Maximize  w  *  ub  +  vb 


(32) 


_  *11  *  *21  1 
subject  to  uA  +  vA  <  c 


22  2 
v  A  <  c 


where  decomposition  is  done  with  respect  to  the  coefficient  vectors  of  the 
variables  v  ,  the  coefficient  vectors  of  the  component  u  are  retained 
unmodified : 


(33) 


Maximize  w  =  ub^  +  £  ARpR  +  1  wk^k 
subject  to  uA11  +  [  A^  +  £  URTR  <  c1 

Ik  -  i 


xk  - 0 


bk  > 0 


k  1c,  2  k 

=  no  for  a  vertex  n  of  the  convex  polyhedron 


where  =  tt^A^  and  p 

22  2  k  21  k  2  k 

determined  by  uA  <  c  ;  TR  =  a  A  and  >R  =  a  b  for  an  extreme  ray  n  of 

22  2 

the  convex  polyhedron,  hA  <  c  .  If  we  now  take  the  dual  of  (33)  assigning 
dual  multipliers  x^  to  the  first  n^  inequalities  and  <’  to  the  last  equation. 
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we  obtain  the  dual  problem: 


(34) 


or  equivalently, 


Minimize  z  **  c^x  +  0 


subject  to  A^x  -  b^" 


V  +  0  >  pk 

-  Yk 


V 


x  >  0 


Minimize  z  «  c^x  +  P 


k  -  1,  •••,  t 
k  "  1,  »  .  .  ,  s 


.  A11 

to  Ax 

=  b1 

,  kA21N 

(o  A  )x 

>  (oV) 

k  -  1 ,  ....  s 

(iA^)x  +  0 

>  (A2) 

k  -  1 . t 

x  >  0 


which  corresponds  to  (22).  Note  that  the  feasibility  constraints  (22. b)  correspond 

22  2 

to  the  extreme  rays  of  the  polyhedron  ttA  <  c  ,  whereas  the  optimality  constraint 

22  2 

(22. c)  correspond  to  extreme  points  of  tiA  <  c  The  constraints  generated  in 
Step  2'  of  the  cutting  plane  algorithm  correspond  to  columns  of  (33)  generated 
during  the  Phase  I  of  this  partial  decomposition  procedure.  Thus,  the  algorithm 
which  we  developed  here  can  be  interpreted  as  a  dual  method  of  the  Dantzig  -  Wolfe 
decomposition  algorithm. 

On  the  other  hand,  if  we  consider  the  L-shaped  linear  program  in  the 


equivalent  form 


Minimize 


c^x  4-  Q(x) 


(35)  subject  to  A^x  - 

x  c  K2 
x  >  0  , 

then  our  algorithm  can  be  interpreted  as  a  cutting  plane  algorithm  [A]  ,  [16]  .  I 

21  22  2 
A  and  A  have  a  finite  number  of  rows,  K  is  a  polyhedral  set  and  Q  a 

convex  polyhedral  function.  The  methods  of  [A],  [16]  can  be  used  to  establish 

the  convergence  of  our  algorithm  in  the  case  where  the  number  of  rows  are 

infinite;  alternatively,  the  results  in  [23]  can  be  used  to  establish  convergence 

using  the  interpretation  of  our  algorithm  as  the  dual  of  a  decomposition  procedure 

Thi:  is  simply  a  reflection  of  the  fact  that  the  cutting  hyperplane  methods 

of  Cheney  and  Goldstein  [A],  Goldstein  [13],  and  Kelley  [16]  on  one  hand  and  the 

decomposition  methods  of  Dantzig  and  Wolfe  [8],  the  algorithm  associated  with 

Wolfe's  Generalized  Program  [5],  [23],  and  in  particular  Dantzig's  convex 

programming  algorithm  [5]  on  the  other  hand  are  simply  dual  methods  to  one  another 
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4.  Optimal  Control  with  State  Constraints 

A  rather  standard  optimal  control  problem  is 


Maximize  q^OF) 

subject  to  ^  -  B(t)q(t)  +  C(t)u(t) 

(36)  q  (0)  -  q° 

T 

q (T)  e  L  =  {q  =  (q0>  ....  qn> | q±  -  qA  i=l,  n) 

q ( t)  e  Q(t)CRn 
u(t)  e  U(t) 

where  U(t)  and  Q(t)  are  closed  convex  polyhedral  sets.  We  will  consider  the 
discrete  analogue  of  this  system. 


Maximize 
subject  to 


N 

<0 

i+1  i 
3 _ I3_ 

A 

0  o 


q 


N 

q 


i 

q 


i 

u 


=  q 

E  L 

C  Ql 


+ 


i 

u 


i  *  0 ,  . . . ,  N-l  where  A 
Since  q*+^  -  [I  +  AB*]q* 


T  i 

=  —  ,  q  =  q(iA)  and  similarly  for  the 
+  AC*u*  ,  we  can  now  solve  for  each  q^ 


other  functions, 
inductively  in 


The  case  where  U(t)  are  not  polyhedral  leads  to  algorithms  which  converge  but 
are  not  finite.  Problems  for  which  U(t)  is  not  polyhedral  is  treated  in  [2'J|. 
Problems  for  which  Q(t)  is  not  polyhedral  can  be  treated  by  analogous  devices. 


o  .  .  ,  0  k-1 

terms  of  the  initial  state  q  and  the  control  sequence  u  ,  ....  u 


q1  -  [I  +  AB° ] q°  +  AC°u° 

q2  -  [I  +  AB1 ] { [ I  +  AB° ]q°  +  AC°u°}  +  AC1u1 


and  in  general 


-  J  ^ nk [ I  +  Afr1]  J  q°  +  [I  +  ABk]  ...  [1  +  AB1]AC°u°  + 


k-1  ,  . „k  k 

4-  [I  +  AB  ] AC  u  +  AC  u 


Y [ j ,k]  -  [I  +  ABk_1)[I  +  ABk'2]  +  ...  +  [I  +  ABJ]  , 


k  R  ^ 

for  j  <k  ,  y [ J . j  1  *  I  and  AC  =  E  .  Then,  we  have 


(38)  qk  -  Y[0,k]q°  +  j  Y[j  +  l,k]EJuj  . 

j-o 

Since  Q*  and  U*  are  closed  convex  polyhedral  sets,  we  can  formulate  the 
constraints  on  in  the  form  F^u^  >  j  =0,  ...»  N-l  ,  and  those  on 

the  state  variables  q  as  C.^^q^  s  .  So  now  we  have 


Maximize 

N-l 

subject  to  U„q^  -  £  Y[j+l,N]E^u^  =  Y[o,N]q 


f(1)..j  ,  f<j> 


i  =  o . N-l 


and  the  additional  constraints 


G^q-^  >  j  *  0 . N  -  1  ,  where  qT  =  (0,  q^,  ...,  q^) 
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If  we  let  u  -  [u° . uN-1]  and  A  =  [-Y(1,N)E°  ,  -Y(2,N)E1 . -Y(N,N)Fk) 

o  T 

and  b  ■*  Y[0,k]q  -  q  ,  we  have 


Maximize 

N 

qo 

(39) 

subject  to 

uoqo  +  Au 

=  b 

F^u^ 

:  f<J) 

j  =  0,  ....  N  -  1 

G^  q^ 

:«(J) 

1=0,  ....  N  -  1 

where  is  given  by  (38).  The  approach  for  handling  the  constraints 

>  f^  by  generalized  linear  programming  has  been  described  in  [6]  and 
[23].  Thus  for  simplicity,  we  limit  ourselves  to  a  discussion  of  the  constraints 
G^q^  >  g^  and  assume  that  the  constraints  F^^u^  >  f^  on  the  u^'s  simply 
reduce  to  the  requirements  that  they  are  all  nonnegative.  We  may  now  simplify  (39) 
to  read: 


Maximize 

(AO)  subject  to 


It  is  this  problem  which  we  interpret  as  an  L-shaped  program.  The  correspondence 
is  A  ~  ,  u  ~  :c  ,  finally  the  slack  variables  of  the  implicit  constraints 

on  u  ,  G^^q^  >  g^^  ,  correspond  to  y  .  In  this  case,  c^  *=  0  so  that 

second  stage  feasibility  is  the  only  requirement.  Frequently  from  the  physical 
nature  of  the  problem  it  is  clear  that  "usually"  the  state  constraints  will  not  be 


N 

q() 


U0q0  +  Au  =  b 

C<jV  >  .(1) 


0  . 
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violated,  and,  of  course,  the  value  of  the  slack  variables  are  of  no  particular 
use  so  that  the  representation  of  (40)  as  an  L-shaped  program  seems  particularly 
appropriate. 

To  generate  the  cut  we  simply  evaluate  it^  by  «  1  for 
[g^^  -  0  and  nj  =  0  otherwise.  The  cut  is  equal  to 

V  n-’c^q^  ■*  ^  ^  which  is  the  sum  of  the  infeasible  equations.  All  that 

remains  is  to  express  these  in  terms  of  the  u's.  In  other  words,  we  wish  to 
evaluate 


(41) 


l  A(k) 


y  Y[j+l,k]EJuj 
j-U 


and  the  constant  term 


I"k[g(k)  -  Y (0,k  )q°]  . 

This  will  give  a  new  constraint  A  , .  . u  b  ,,  which  must  be  satisfied  by  u  , 

where  An+^  *  denotes  the  n+1^  row  of  A  ,  the  special  structure  of  (41),  in 

particular,  of  the  Y(j,k]  makes  possible  many  simplifications  in  determining 

^n+1*  and  ^n+1  and,  in  particular,  the  relevant  quantifies  would  be  accumulated 

as  one  determines  ,  rather  than  determining  ti^  and  then  going  back  to 

calculate  A  , .  and  b  ,  .  In  addition,  if  the  state  space  constraints  are 
n+1*  n+1 

"loose",  not  many  of  the  equations  would  be  violated. 

This  application  is  an  example  of  an  important  subclass  of  L-shaped  programs 
which  could  be  called  I-shaped  programs.  These  are  L-shaped  programs  in  which  the 
components  of  the  v  vector  are  simply  slack  variables. 

The  'nteger  programming  algorithm  of  Gomory  [14]  can  be  considered  as  another 

21  2 

example  of  an  I-shaped  program  where  A  x  +  Iv  =  b  or  equivalently 
1 1  2 

A  x  <  b  represents  the  infinite  number  of  constraints  which  can  be  added  to 
eliminate  noninteger  extreme  points  but  do  not  eliminate  anv  feasible  integer  points. 
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3.  Stochastic  Programs  with  Recourse 

A  stochastic  program  with  recourse  (random  right-hand  sides)  also  known  as 
two  stage  linear  programs  under  uncertainty  [7]  reads 

(42)  Minimize  z  =  cx  +  Er(Min  qv) 

subject  to  Ax  =  b 

(42. a)  Tx  +  Wy  >  £  f,  on  (c.z.F) 

x  >  0  ,  y  >  0 

The  interpretation  to  be  given  to  this  problem  as  well  as  the  definition  of  the 
symbols  can  be  found  in  [25]  or  various  other  papers  in  this  area,  see  e.g.  [7], 

[  1  5]  -  Problem  (42)  is  easily  recognized  to  be  an  I. -shaped  program  with  possibly 
an  infinite  number  of  constraints  (42. b)  and  an  infinite  number  of  v  variables. 
We  denote  by  the  support  of  the  random  variable  :  ,  i.e.,  the  smallest  closed 

subset  of  Rm  of  measure  one. 

V'e  shall  assume  that  -  has  at  least  upper  bound  a  ,  i.e.,  sucli  that 

ac=  and  for  all  i  ,  f  <  for  all  f,t"  .  If  this  model  is  viewed  as  the 

representation  of  a  physical  decision  process,  the  assumption  that  for  each  i 

there  exist  a.  such  that  f. .  ■  <».  seems  to  be  very  natural.  The  additional 
l  l-i 

assumption  that  ar.~  is  somewhat  more  restrictive.  However,  this  would 
certainly  be  the  case  if  the  components  of  >.  where  independent  random  variables 
and  each  t  has  compact  (or  bounded  above)  support.  Extensions  and  a  more 
complete  discussion  of  these  questions  can  be  found  in  "Finding  a  Feasible  Solution 
to  Stochastic  Program  with  Fixed  Recourse"  (24). 

From  a  mathematical  viewpoint  the  assumption  that  for  each  i  there  exists 


i  . 
l 


is  not  so  appealing  but  if  such  an  upper  bound  does  not  exist  then  determininr. 
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if  problem  (42)  is  feasible  has  to  be  dealt  with  differently,  as  can  be  seen  from 
the  following  proposition. 


(43) 


Proposition : 


Suppose  for  some  i  there  is  no  number  such  that 

£  <  and  £eE  .  Then  (42)  is  feasible  only  if  the  lir.oality 


riin 


m 


space  of  pos  (W,  -I)  contains  R^  ,  where  R^  is  the  i 
component  of  the  Cartesian  product  R1 
real  line). 


th 


2iR  (Rj  denotes  the 


Proof : 

If  R^C  pos  (W,  -I)  then  the  equation  (W,  -1)^  (^)  **  4^  Is  solvable  for  all 
£  with  y  and  s  nonnegative.  Otherwise,  for  some  4^  the  above  equation 
is  not  solvable.  Since  £^  has  no  upper  bound,  for  any  x  there  exists  £ 

In  5  (determining  £^)  such  that  the  system 


V  >  V 

y  >  0 

is  inconsistent.  This  implies  that  for  no  x  the  recourse  (or  second  stage) 
problem  is  feasible  for  all  £  in  E  ,  thus  the  set  of  feasible  solutions  to 
(42)  is  empty,  i.e.,  (42)  is  infeasible. 

If  the  £j's  are  independent  and  for  certain  i's  ,  £^  has  no  greatest 
upper  bound,  we  can  use  the  above  proposition  to  determine  if  (42)  is  infeasible. 

If  the  criterion  is  satisfied  we  can  ignore  those  equations  whenever  we  verify  if 
a  given  x  is  feasible  or  not.  In  (28]  the  problem  of  characterizing  and  computing 
(which  can  be  easily  done)  the  lineality  space  of  pos(W,  -1)  has  been  dealt  with 


in  detail. 
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In  the  algorithm  to  be  described  below,  we  shall  assume  that  for  each  i  , 
exists  and  aeE  .  Propsoition  (16)  in  Section  2.B  "A  Feasibility  Test",  of 
[26]  allows  us  to  derive  constraints  (equations  (17)  of  [26]  of  the  form) 

(44)  (cT)x  ^  oa 

which  in  view  of  proposition  (16)  of  [26],  plays  the  same  role  than  the  feasibility 
constraints  (7)  plays  in  the  L-shaped  linear  program.  Moreover,  it  has  been  shown 
that  the  feasibility  region  for  the  decision  variables  x  determined  by  the  induced 
constraints  [26,  p . 9 2 ]  can  be  represented  by  a  finite  number  of  linear  constraints. 
(Prosposition  12  of  [27].)  In  Section  2.D  "Some  Geometric  Characterizations,"  we 
have  shown  the  relation  between  the  feasibility  constraints  (7)  that  we  introduce 
and  the  supports  of  the  cone  pos  (W,  -I)  .  In  [27]  the  accent  had  been  placed 
on  deriving  an  expression  in  terms  of  a  minimal  number  of  support  (determined  by 
the  rows  of  the  polar  matrix  [27])  of  pos  (W,  -I)  ,  rather  than  an  arbitrary 
finite  collection  of  supports.  As  can  be  seen  from  proposition  (24)  supports  of 
maximum  dimension  corresponds  to  obtaining  a  particular  solution  to  the  lineai 
program. 


Minimize  ev 

+  -  i 

subject  to  Wy  +  Iv  -  Iv  =  a  -  Tx 

y  >  0  ,  v+  >  0  ,  v  >0. 

These  observations  allow  us  to  construct  an  algorithm  which  will  find  feasible 
solutions  to  (42)  in  a  finite  number  of  steps,  i.e.,  by  requiring  that  x  satisfies 
the  constraints  (42. a)  to  which  we  add  a  finite  number  of  constraints  of  the 
form  (44);  each  one  being  generated  by  solving  om  linear  program  rather  than 
verifying  it  for  a  particular  x  and  for  all  there  exist  a  feasible  y  , 
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i.e.  ,  y  >  0  ,  such  that  Wy  >  E,  -  Tx  . 

We  now  outline  a  general  algorithm  for  solving  problem  (42);  general,  in 
the  sense  that  we  make  no  assumption  on  the  structure  of  the  matrices  (in 
particular  W)  or  on  the  form  of  the  distribution  of  the  random  variable  £  , 
except  that  5  has  a  greatest  upper  bound.  (See  proposition  (43)  if  this  is 
not  the  case.)  We  ignore  the  special  cases  of  infeasibility  and  unboundedness 
which  are  to  be  handled  as  before. 

Step  1: 

Solve  the  linear  program 


(45) 

Minimize  cx  -*■ 

9 

(45 . a) 

subject  to  Ax 

(45. b) 

(okT)x  > 

Q 

7? 

P 

k  -  1 . s 

(45. c) 

(TrkT)x  + 

9  >  pk 

k  »  1 ,  .  .  .  ,  t 

x  >  0  . 

Initially,  s  and  t  are  zero.  If  no  constraints  of  the  form  (45. c)  are  present 

l  l 

9  is  set  equal  to  -  00  and  is  ignored  in  the  computation.  Let  x  ,  6  be  an 
optimal  solution  of  (45). 

Step  2: 

Solve  the  linear  program  to  find 


w  ^  -  Minimum  ev+ 


(46) 


+  -  £ 

subject  to  Wy  +  Iv  -  Iv  -  a  -  Tx 
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11  £ 

If  w  =  0  ,  go  to  Step  3.  If  w  4  0  ,  the  optimal  solution  o  to  (46)  is  used 

to  generate  a  cut  of  the  form  (45. b). 

Step  3: 

For  all  in  E  ,  solve  the  linear  program. 


w  =  Minimum  gy 


(47) 


subject  to  Wv  -  Is  =  >.  -  Tx 


y  >  0  ,  s  >  0  . 

I 

Each  £  determines  an  optimal  it  ,  say,  it  (£,)  .  We  then  compute 

w^(x2)  -  Er  { TT  2  (O  (£;  -  Tx2)}  ,  tt2  =  Ef.{TT2(0}  and  p2  =  Ef  (  tt  2  (f,)f, }  .  If 

2  £  £  £  £ 
w  (x  )  <  0  ,  we  terminate  (Optimality  Criterion).  If  not,  we  use  tt  ,  p  to 

generate  a  new  constraint  of  the  form  (45. c)  which  we  now  add  to  our  problem  (45) 

and  return  to  Step  1. 

We  should  also  point  out  that  in  following  this  procedure,  it  is  possible 
to  generate  an  infinite  number  of  constraints  of  the  form  (45. c;.  Nevertheless,  a 
result  of  K.  Murty  [17]  allows  us  to  keep  m  (T  is  m  x  n)  or  less  constraints 
of  the  form  (45. b)  and  (45. c)  at  each  cycle,  i.e.,  the  constraints  with  nonzero 
slack  can  be  removed. 

We  have  separated  Step  2  of  the  paraphrase  (in  Section  21, c)  of  the  cutting 
plane  algorithm  in  two  parts.  The  reason  being  that  in  order  to  generate  the 
feasibility  cuts,  we  need  only  consider  the  upper  bound  of  E  (not  all  elements 

of  E)  whereas  we  need  complete  information  related  to  the  probabilitv  space 

£  £  I  . 

(- ,  F,  F)  in  order  to  compute  tt  and  p  (even  when  has  finite 
cardinality  the  labor  so  saved  should  be  considerable).  Moreover,  if  has 

infinite  cardinality,  it  is  difficult  to  perform  Step  3  unless  the  structure  of 

f  f 

W  is  such  that  it  is  possible  to  find  a  closed  form  expression  for  n  and  p  , 


and 
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e.g.,  see  [20]  and  [25].  The  remaining  part  of  this  section  is  denoted  to  suggest 
a  method  to  circumvent  this  problem.  We  start  by  describing  a  variant  of  the 
above  algorithm. 

If  £,  is  an  absolutely  continuous  random  variable,  we  can  modify  the 
algorithm  as  follows: 

Step  1: 

Solve  the  linear  program 


(48) 


Minimize 


subject  to 


[c  -  IT 


l-l 


A11  k 

A  x  -  b 


(0* T)x  > 


x  > 


T  ]  x 

(ok  a) 
0  . 


k  *  1 ,  ....  s 


H-l  -£ 

Initially  s  *  0  and  it  =  0  .  Let  x  be  an  optimal  solution  to  (48).  Find 

Min  U/ ( X )  =c[(l-A)x^  *  +  Ax^]  +  Q[(l-A)x^  *  +  Ax^]  where  x^  *  was  our  previous 

0  £-1 
solution  (x  *  0)  which  for  £  >  1  was  used  to  determine  tt  and  the 

t 

function  Q(x)  is  as  defined  in  [26,  Equation  (21)].  Say,  iJ>(A  )  <  4<(A)  for 

£  £-1 
A  e  [0,1].  If  A  =  0  .  We  terminate  with  optimal  solution  x 

(Optimality  Criterion).  Let  x^  =  (1  -  A^)x^  *  +  A^x^  . 


Step  2: 
As  above. 


Step  3: 

l 

As  above,  determines  ti  and  we  npu  return  to  Step  1. 

The  convergence  of  this  algoritL-  can  be  easily  verified  if  we  observe  that 
from  proposition  (29)  and  corollary  (28)  in  [26]  it  follows  that  if  f.  is  a 

continuous  random  variable  then  Q(x)  is  a  differentiable  function  with  gradient 

£  } 

-ti  T  at  x  Thus,  the  above  algorithm  can  be  viewed  as  a  variant  nl  the 
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Frank-Wolfe  [10]  algorithm  for  finding  the  minimum  of  a  convex  differentiable 

function  on  a  convex  polyhedral  set.  Since  their  procedure  requires  twice- 

2 

differentiability  and  in  general  Q(x)  will  not  be  of  class  C  their 
convergence  proof  does  not  strictly  apply.  The  necessary  modifications  can  be 
found  in  [25,  proposition  (37),  (41)  and  (43)]. 

This  last  algorithm  as  well  as  the  first  one  we  suggested,  to  solve 
Stochastic  Programs  with  Recourse  relies  on  the  possibility  of  performing  Step  3  . 
If  5  does  not  have  finite  cardinality  this  seems  to  be  nearly  impossible. 
However,  one  could  exploit  a  suggestion  of  Dantzig  and  Madanskv  [7]  which  consists 

in  sampling  the  distribution  of  £,  and  solve  Step  3  for  some  finite  sample.  This 

i  £ 

would  naturally  result  in  approximated  values  for  it  and  p  .  As  has  already 
been  pointed  out  in  the  introduction,  this  approach  would  generate  a  constraint 
of  tvpe  (45. c)  which  would  not  necessarily  be  a  support  of  the  function  Q(x)  , 
and  could  bery  well  eliminate  the  optimal  solution  (42)  from  the  set  of  feasible 
solutions  to  (45).  This  inconvenience  has  been  completely  eliminated  if  we  follow 
the  second  procedure  since  all  the  constraints  present  in  (48)  never  involve  any 
approximation  process. 

We  are  however  still  left  with  two  problems.  First  to  solve  Step  1  for  a 

large  (possibly  very  large)  number  of  values  of  f.  in  -  .  Second,  the 

'  i  ? 

resulting  tt  will  not  in  general,  determine  the  gradient  of  Q(x)  at  x  and 

thus  the  convergence  properties  of  the  algorithm  are  changed.  This  second  problem 
will  be  che  object  of  another  paper  in  which  various  sampling  techniques  are 

examined  and  the  convergences  properties  of  the  algorithm  are  established.  We 

'  i  t 

now  show  how  to  obtain  the  approximate  value  tt  for  tt  from  a  specific  sample. 
1  N 

.  et  t,  ,  .  .  .  ,  f,  be  a  sample  of  size  N  obtained  from  the  distribution  of  f  . 

Our  purpose  is  to  solve  (in  Step  3)  the  N  linear  programs  of  the  form: 
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Minimize  qy 


(49) 


subject  to  Wy  -  Is 


k  -  1 . N 


y  >  0  s  >  0 

Since  we  are  performing  Step  3,  x  c  K  ,  and  since  ^  t  ■  it  follows  that  (49) 

i  It 

is  feasible  for  all  £  .  Moreover  if  (49)  is  unbounded  for  some  f,  it  is 
unbounded  for  all  £  e  5  thus  (42)  is  unbounded  (see  Proposition  13).  Let  us 

p 

assume  (49)  is  solvable.  Let  n  (£)  denote  the  optimal  simplex  multipliers 
associated  with  solving  (49)  for  a  particular  £  .  Then 


(50) 


'l 

71 


H  1 

■  I  h 

k-1  w 


1  uk) 


In  the  appendix  of  [26],  we  have  reviewed  the  properties  of  the  function 
Q(t)  *=  (Min  qy|Wy  ■  t,  y  >  0}  in  particular  we  shall  use  the  fact  that:  If 

U(t)}  is  the  mapping  determining  the  optimal  simple*  multipliers,  then  there 

exists  a  function  7T(t)  in  { tt ( t ) }  piecewise  constant  on  pos  W.  In  particular,  if 
is  an  optimal  basis  corresponding  to  a  particular  value  of  t  ,  then  is 

also  an  optimal  basis  for  all  t  e  pos  ^  .  Let  be  the  subvector  of  q 

corresponding  to  W^11  ,  then  tt ( t )  ■  q^W^  1  determines  an  optimal  vector 
of  multipliers  for  all  t  e  pos  .  Applying  this  result  to  (49)  it  follows 

that  w  (£)  is  also  optimal  for  all  £^  such  that  (£^  -  Tx*)  e  pos  W(£)  where 
W(£)  is  the  optimal  basis  obtained  from  solving  (49)  for  some  fixed  £  .  To 
determine  if  (£^  -  Tx  )c  for  W(£)  it  is  sufficient  to  verify  if 


W(0-1  (£j  -  Tx4)  >  0  . 


-1 


This  can  be  easily  done  since  W(£)  is  available  from  the  final 


We  now  give  an  algorithmic  procedure  to  find  ttX  ,  as  defined 


optimal  tableau, 
by  (50). 


k. 
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Step  a: 

Select  an  unbiased  sample  of  size  N  from  the  distribution  of  (,  ,  sav 

IN  k  k  l  1 

S  ,  £,  .  Compute  r,  =  £  -  Tx  k  =  1 ,  ....  N  .  By  { }  we  denote  the 

set  of  available  and  set  L  =  N  . 

Step  b: 

1  k 

Select  some  c,  in  { £,J }  and  set  c,  =  c,  (initially  k  =  1)  and  solve 
the  linear  program: 


Minimize  qy 


subject  to  Wy  -  Is  =  r. 


y  >  0  . 


k  k 

Let  )  be  the  optimal  simplex  multipliers  and  W(C  )  be  the  corresponding 

optimal  basis. 

Step  c: 

Let  n(k)  be  the  number  of  vectors  r,^  in  the  set  { ^ }  such  that 


(51)  WUV1  >  0  . 


Set 

L  =  L  -  n(k) 

and 

if  L  >  0 

return  to  Step  b  with 

k  =  k  +  1 

and  delete 

from 

the  set  {c^} 

and 

those 

which  satisfied  (51). 

If  L  =  0 

terminate 

with 


i 

7? 


\ 

N 


I  n(k)n (rk)  . 
k 


of 


In  returning  to  Step  b  it  is  suggested  to  select  c,  (in  the  remaining  set 
{^}  such  that  r,  fails  to  satisfy  (51)  only  in  a  minimum  number  of 
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components  (if  possible  one).  Thus  the  previous  basis  would  be  the  optimal  basis 
k 

for  the  new  c,  up  to  very  few  dual  simplex  step. 

A  few  experiments  have  been  made  on  an  IBM  7094  (with  a  not  nearly  optimal 

code).  We  have  selected  N  =  3000  and  5000,  and  10  <  m  <  40  (m  is  the  numbers 

'  £ 

of  rows  of  W  ).  In  each  case  the  computation  of  n  took  never  more  than  twice 
the  time  required  to  solve  one  linear  program  of  the  same  size.  In  the  same  vein, 
a  number  of  experiments  have  been  conducted  by  Ballinfty  and  Prekopa  for  random 
linear  programs.  In  their  manuscript  [1]  they  show  that  numerous  "tricks"  can  be 
performed  to  improve  sampling  procedures. 
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